function [Q,errorMes,model,setupFilter,outFilter] = objectFuncQML(paramsInputValues,setupFilter)
Q         = NaN;
model     = [];
outFilter = [];
% Check if the lower and upper bounds are violated
if isfield(setupFilter,'transformParamsOn') == 1
    paramsInputValues = unTransformParams(paramsInputValues,...
        setupFilter.lowerBoundsValues,setupFilter.upperBoundsValues);
else
    if any(paramsInputValues < setupFilter.lowerBoundsValues) ...
            || any(paramsInputValues > setupFilter.upperBoundsValues)
        errorMes = 1;
        return;
    end
end

% We turn paramsInputValues into a struct 
for i=1:length(setupFilter.selectParams)
    name = setupFilter.selectParams(i,1);
    paramsInput.(name{1}) = paramsInputValues(i,1);
end

% Accounting for calibrated coefficients, we unold paramsInput by to get params 
params = paramsInput;
% The setdiff implies that we do not overwrite elements in params based on calibrateparams
namesCalibrate = setdiff(fieldnames(setupFilter.calibrateParams),setupFilter.selectParams);
for i=1:length(namesCalibrate)
   name = namesCalibrate(i);  
   params.(name{1}) = setupFilter.calibrateParams.(name{1}); 
end

% Setting the size of the measurement errors
params.stdMhours =  params.stdMhours*params.ScalingMacro;
params.stdMwage  =  params.stdMwage*params.ScalingMacro;
params.stdMdc    =  params.stdMdc*params.ScalingMacro;
params.stdMinf   =  params.stdMinf*params.ScalingMacro/2;

params.stdMr1    =  params.stdMr1*params.ScalingYield;
params.stdMr4    =  params.stdMr4*params.ScalingYield;
params.stdMr12   =  params.stdMr12*params.ScalingYield;
params.stdMr20   =  params.stdMr20*params.ScalingYield;
params.stdMr28   =  params.stdMr28*params.ScalingYield;
params.stdMr40   =  params.stdMr40*params.ScalingYield;
params.stdSurvey =  params.stdSurvey*params.ScalingYield;

if ~isfield(params,'ZI')
    params.ZI = (1-params.THETA+params.ETA*params.THETA)*(params.ETA-1)*params.XI/...
        ((1-params.XI)*(1-params.THETA)*(1-params.XI*params.BETTA*params.MUZss^(1-params.CHI*(1-params.CHI0)-params.CHI0)));
end

if ~isfield(params,'ALFA') && (isfield(params,'U0d') || isfield(params,'U0'))
    params.ALFA = getALFA(params);
    if abs(params.ALFA) > 500
        errorMes = 1;
        return;
    end
elseif ~isfield(params,'U0d') && isfield(params,'U0')
    params.U0d = getU0d(params);
elseif isfield(params,'U0d') && ~isfield(params,'U0')
    params.U0 = getU0d(params);
end

%% Solve the model
% Starting values from a perturbation approximation
setupProj   = setupFilter.setupProj;
appOrderPer = 2;
[modelPer,errorMes] = perturbationDSGEmodel(params,appOrderPer,...
    setupProj.setupModel,setupFilter.MEXon*0);
if errorMes == 1
    disp('DSGE model: No unique and stable solution')
    Q = NaN;
    errorMes = 1;
    return
end
modelPer.inclSurveys = setupFilter.inclSurveys;

% Version 1: from a log-linear approximation
varX          = reshape((eye(modelPer.nx^2)-kron(modelPer.hx,modelPer.hx))\reshape((modelPer.eta*modelPer.eta'),modelPer.nx^2,1),modelPer.nx,modelPer.nx);
diagStdX      = max(sqrt(diag(varX)),0.001);
Smol_grid_iso = setupProj.Smol_grid_iso;
%xGrid         = Smol_grid_iso  .*repmat((setupProj.gridFactor_x*diagStdX'),size(Smol_grid_iso,1),1);
meanXPer      = (eye(modelPer.nx)-modelPer.hx)\(0.5*modelPer.hss+0.5*reshape(modelPer.hxx,modelPer.nx,modelPer.nx^2)*varX(:));
xGrid         = Smol_grid_iso.*repmat((setupProj.gridFactor_x*diagStdX'),size(Smol_grid_iso,1),1);
if setupFilter.setupProj.endoMeanGrid == 0
    xGrid = xGrid + repmat(meanXPer',size(Smol_grid_iso,1),1);
end

setupProj.appOrderPer = appOrderPer;
setupProj.modelPer    = modelPer;
%setupProj.numCPUs     = numCPUs;
setupProj.gSS         = modelPer.g0;
setupProj.hSS         = modelPer.h0;
setupProj.hx          = modelPer.hx;
setupProj.eta         = modelPer.eta;
setupProj.params      = modelPer.params;
setupProj.nx1         = modelPer.nx1;
setupProj.nx          = modelPer.nx;
setupProj.ny          = modelPer.ny;
setupProj.mx          = modelPer.mx;
setupProj.myx         = modelPer.myx;
setupProj.labelx      = modelPer.labelx;
setupProj.labely      = modelPer.labely;
setupProj.scaleX_in_g = diagStdX;
setupProj.xGrid       = xGrid;

fileStart = ['startValues_appOrder',num2str(setupProj.appOrder)];
[theta0,vecTheta0,setupProj] = loadThetaStart(fileStart,setupProj);
setupProj.nDim1Theta= size(theta0,1);
setupProj.nDim2Theta= size(theta0,2);

%setupProj.objectFuncNum = 1; %for main model
%[Q0,resEuler0] = NLSobjectFunc(vecTheta0,setupProj);
%Q0'*Q0

setupProj.objectFuncNum = 1; %main model
modelNoYields = projectionStep(vecTheta0,setupProj);
if modelNoYields.exitflag <= 0
    disp('DSGE model: No Projection Solution')
    Q = NaN;
    errorMes = 1;
    return
end

if setupFilter.setupProj.endoMeanGrid == 1
    setupFilter.setupProj.xGrid = modelNoYields.xGrid; %update of grid
end
modelNoYields = ProjLoadingsExtraControls(modelNoYields);
setupProj.objectFuncNum = 2; %for bond prices
setupProj.maxMat = max(setupFilter.maturities); 
%model = projectionStepBondPrices(modelNoYields,setupProj);
model = projectionStepBondPricesClosedForm(modelNoYields,setupProj);

if setupFilter.inclSurveys == 1
   % Solve for expected short rates
   var0   = model.by0(1,1);
   varx   = model.byx(1,:);
   varxx  = model.byxx(1,:,:);
   varxxx = model.byxxx(1,:,:,:);
   %[r0,rx,rxx,rxxx]= projectionStepExpectedControl(var0,varx,varxx,varxxx,...
   %    max(setupFilter.matConShortRate),model,setupProj);
   [r0,rx,rxx,rxxx]= projectionStepExpectedControlClosedForm(var0,varx,varxx,varxxx,...
       max(setupFilter.matConShortRate),model,setupProj);
   model.matConShortRate = setupFilter.matConShortRate;
   model.r0     = r0;
   model.rx     = rx;
   model.Rxx    = reshape(rxx,length(setupFilter.matConShortRate),model.nx^2);
   model.Rxxx   = reshape(rxxx,length(setupFilter.matConShortRate),model.nx^3);
end
model.maturities  = setupFilter.maturities;
model.inclSurveys = setupFilter.inclSurveys;
model.pruningOn   = setupFilter.pruningOn;

%Updating setupProj in setupFilter
setupFilter.setupProj = setupProj;     
%% Running the filter
% Setting the co-variance matrix for the measurement errors
if setupFilter.inclSurveys == 1
    dimy  = setupFilter.numMacro+length(model.maturities)+length(model.matConShortRate); 
else
    dimy  = setupFilter.numMacro+length(model.maturities);
end
Rw        = zeros(dimy,dimy);
Rw(1,1)   =  params.stdMhours^2;
Rw(2,2)   =  params.stdMwage^2;
Rw(3,3)   =  params.stdMdc^2;
Rw(4,4)   =  params.stdMinf^2;
Rw(5,5)   =  params.stdMr1^2;
Rw(6,6)   =  params.stdMr4^2;
Rw(7,7)   =  params.stdMr12^2;
Rw(8,8)   =  params.stdMr20^2;
Rw(9,9)   =  params.stdMr28^2;
Rw(10,10) =  params.stdMr40^2;
if setupFilter.inclSurveys == 1
    Rw(11:dimy,11:dimy) = params.stdSurvey^2*eye(dimy-10);
end
model.Rw= Rw;
if setupFilter.pruningOn == 1
    % Initializing the filter
    % We use the pruned perturbation approximation
    vec_xfxf  = reshape(varX,model.nx^2,1);
    if model.appOrder == 1
        xf0       = (eye(model.nx)-model.hx)\model.h0;
        xAll0     = xf0;
        PxAll0    = ones(model.nx+0*model.nx1,model.nx+0*model.nx1)*1D-5;
        PxAll0(1:model.nx,1:model.nx) = reshape(vec_xfxf,model.nx,model.nx);
    elseif model.appOrder == 2
        xf0       = (eye(model.nx)-model.hx)\model.h0;
        xs0       = (eye(model.nx)-model.hx)\(model.Hxx*(vec_xfxf+reshape(xf0*xf0',model.nx^2,1)));
        xAll0     = [xf0;xs0(1:model.nx1,1)];
        PxAll0    = ones(model.nx+1*model.nx1,model.nx+1*model.nx1)*1D-5;
        PxAll0(1:model.nx,1:model.nx) = reshape(vec_xfxf,model.nx,model.nx);
    elseif model.appOrder == 3
        xf0       = (eye(model.nx)-model.hx)\model.h0;
        xs0       = (eye(model.nx)-model.hx)\(model.Hxx*(vec_xfxf+reshape(xf0*xf0',model.nx^2,1)));
        xrd0      = xf0*0;
        xAll0     = [xf0;xs0(1:model.nx1,1);xrd0(1:model.nx1,1)];
        PxAll0    = ones(model.nx+2*model.nx1,model.nx+2*model.nx1)*1D-5;
        PxAll0(1:model.nx,1:model.nx) = reshape(vec_xfxf,model.nx,model.nx);
    end
    % Running the filter
    [outFilter,errorMes]=CDKF_dsge(model,setupFilter,setupFilter.data(:,1:dimy),xAll0,PxAll0,Rw);
else    
    % Initializing the filter
    x0      = (eye(model.nx)-model.hx)\model.h0;
    P0      = varX;
    
    % Running the filter
    [outFilter,errorMes]=CDKF_dsge(model,setupFilter,setupFilter.data(:,1:dimy),x0,P0,Rw);
end

if errorMes ~= 0
    return;
end
if setupFilter.SMMwithCSreg == 2
    % Adding term premia for nominal yields to model
    model = getTermPremia_loadings(model,setupFilter);
    termPrem = getTermPremia(model,outFilter.xhat);
    % Risk-adjusted CS regression loadings
    % Convert into forwards
    fwdPrem = nan(size(termPrem.TP));
    for i = 1:size(termPrem.TP,1)-1
        fwdPrem(i,:) = (i+1)*termPrem.TP(i+1,:) - i*termPrem.TP(i,:);
    end
    % TPxhat and hence fwdPrem are scaled by 40000 to get risk premia in
    % annualized basis points in getTermPremia.
    % Annualized yields in are not scaled. Hence, we need to undo the scaling of TP
    % to get basis points, but we still need to preserve the annualization i.e. the multiplication
    % by 4.
    numBoots = 0;
    CS_riskAdj = CampbellSchillerRegression_riskAdjust(setupFilter.yieldsAll,...
        termPrem.TP'/40000*4,fwdPrem'/40000*4,setupFilter.CSregress_m,numBoots);
end
%% Unconditional Moments for Shrinkage 
addMerror = setupFilter.addMerror;
if setupFilter.numSim > setupFilter.burnIn
    % Moments computed by simulation
    outSim      = simProjection(model,setupFilter.shocks(:,1:setupFilter.numSim));
    if setupFilter.pruningOn == 1
        %yieldsAllSim   = max(getYields(model,outSim.xAll_cu),0);
        yieldsAllSim   = getYields(model,outSim.xAll_cu);
    else
        yieldsAllSim   = getYields(model,outSim.x_cu);
    end
    % Adding measurement errors
    yieldsAllSim         = 4*yieldsAllSim + addMerror*model.params.stdMr1*setupFilter.noiseYieldsAll(1:size(yieldsAllSim,1),1:size(yieldsAllSim,2));
    yieldsAllSimNoNoise  = 4*yieldsAllSim;
    
    pos_l                = find(strcmp({'$l_t$'},model.labely));
    pos_w                = find(strcmp({'$w_t$'},model.labely));
    pos_c_ba1            = find(strcmp({'$c_{t-1}$'},model.labelx));
    pos_myz              = find(strcmp({'$\mu _{z,t}$'},model.labelx));
    pos_c                = find(strcmp({'$c_t$'},model.labely));
    pos_pai              = find(strcmp({'$\pi_t$'},model.labely));
    ySimAsInData         = zeros(setupFilter.numMacro+length(model.maturities),setupFilter.numSim);
    ySimAsInDataNoNoise  = zeros(setupFilter.numMacro+length(model.maturities),setupFilter.numSim);
    
    ySimAsInData(1,:)    = outSim.y_cu(pos_l,:)-model.gSS(pos_l,1) + addMerror*model.params.stdMhours*setupFilter.noise(1,1:size(outSim.y_cu(pos_l,:),2));
    ySimAsInData(2,:)    = outSim.y_cu(pos_w,:)-model.gSS(pos_w,1) + addMerror*model.params.stdMwage*setupFilter.noise(2,1:size(outSim.y_cu(pos_w,:),2));
    ySimAsInData(3,:)    = 4*(outSim.y_cu(pos_c,:)-(outSim.x_cu(pos_c_ba1,:)+model.gSS(pos_c_ba1,1))+outSim.x_cu(pos_myz,:)+log(model.params.MUZss)) ...
                        + addMerror*model.params.stdMdc*setupFilter.noise(3,1:size(outSim.y_cu(pos_c,:),2));
    ySimAsInData(4,:)    = 4*outSim.y_cu(pos_pai,:) ...
                        + addMerror*model.params.stdMinf*setupFilter.noise(4,1:size(outSim.y_cu(pos_pai,:),2));
    ySimAsInData(setupFilter.numMacro+1:setupFilter.numMacro+length(model.maturities),:) ...
                         = yieldsAllSim(model.maturities,:);
    
    % Without noise added
    ySimAsInDataNoNoise(1,:)      = outSim.y_cu(pos_l,:)-model.gSS(pos_l,1);
    ySimAsInDataNoNoise(2,:)      = outSim.y_cu(pos_w,:)-model.gSS(pos_w,1);
    ySimAsInDataNoNoise(3,:)      = 4*(outSim.y_cu(pos_c,:)-(outSim.x_cu(pos_c_ba1,:)+model.gSS(pos_c_ba1,1))+outSim.x_cu(pos_myz,:)+log(model.params.MUZss));
    ySimAsInDataNoNoise(4,:)      = 4*outSim.y_cu(pos_pai,:);
    ySimAsInDataNoNoise(setupFilter.numMacro+1:setupFilter.numMacro+length(model.maturities),:) ...
                                  = yieldsAllSimNoNoise(model.maturities,:);                 
                     
    burnIn                        = setupFilter.burnIn;
    outFilter.yieldsAllSim        = yieldsAllSim;
    outFilter.yieldsAllSimNoNoise = yieldsAllSimNoNoise;
    outFilter.outSim              = outSim;
    outFilter.ySimAsInData        = ySimAsInData;
    outFilter.ySimAsInDataNoNoise = ySimAsInDataNoNoise;
    if any(any(abs(outSim.y_cu')> 500)) || any(any(isnan(outSim.y_cu')))
        disp('DSGE model: explodes in simulation')
        errorMes = 1;
        return
    else
        if setupFilter.SMMwithCSreg == 2
            [modelSMM,CSmodel] = momentsSMM(ySimAsInData(:,burnIn+1:end)',...
                yieldsAllSim(:,burnIn+1:end)',setupFilter.nyMom,setupFilter.CSregress_m,1);
            % We add the covariances for the risk-adjusted CS moments to the model moments 
            select   = setupFilter.maturites_CS/4;
            T        = size(CS_riskAdj.Ydata,1);
            tmp_yx   = (CS_riskAdj.Ydata(1:T,select)-repmat(mean(CS_riskAdj.Ydata(1:T,select),1),T,1)).*CS_riskAdj.Xdata(1:T,select);
            % For checking
            %tmp_xx   = (CS_riskAdj.Xdata(1:T,select)-repmat(mean(CS_riskAdj.Xdata(1:T,select),1),T,1)).*CS_riskAdj.Xdata(1:T,select);
            %[CS_riskAdj.CSbetta(2,setupFilter.maturites_CS/4) mean(tmp_yx,1)./mean(tmp_xx,1)]
            modelSMM.dataMoments = [modelSMM.dataMoments mean(tmp_yx,1)];
        else
            [modelSMM,CSmodel] = momentsSMM(ySimAsInData(:,burnIn+1:end)',...
                yieldsAllSim(:,burnIn+1:end)',setupFilter.nyMom,setupFilter.CSregress_m,setupFilter.SMMwithCSreg);
        end
        outFilter.modelMoments = modelSMM.dataMoments;
        
        % Extra moments to evaluate the fit of the model
        outFilter.resMoments.CS_loadings  = CSmodel.CSbetta(2,:);
        outFilter.resMoments.meanYields   = mean(yieldsAllSim(:,burnIn+1:end),2);
        outFilter.resMoments.stdYields    = std(yieldsAllSim(:,burnIn+1:end),0,2);
        outFilter.resMoments.meanForEstim = mean(ySimAsInData(1:setupFilter.nyMom,burnIn+1:end),2);
        outFilter.resMoments.stdForEstim  = std(ySimAsInData(1:setupFilter.nyMom,burnIn+1:end),0,2);
    end
elseif setupFilter.tau == 0
    if any(model.heteroShocks ~= 0)
        error('Closed-form solution for GMM moments not implemented with hetero-shocks')
    end
    if model.appOrder == 3
        error('Closed-form solution for GMM moments not implemented at third order')
    end
    % Moments computed in closed form based on the pruned approximation
    outFilter.resMoments   = momentsClosedForm(model,setupFilter);
    scaling                = outFilter.resMoments.scaling;
    if setupFilter.SMMwithCSreg == 2
        % We add the covariances for the risk-adjusted CS moments to the model moments 
        select   = setupFilter.maturites_CS/4;
        T        = size(CS_riskAdj.Ydata,1);
        tmp_yx   = (CS_riskAdj.Ydata(1:T,select)-repmat(mean(CS_riskAdj.Ydata(1:T,select),1),T,1)).*CS_riskAdj.Xdata(1:T,select);
        outFilter.resMoments.modelMoments = [outFilter.resMoments.modelMoments mean(tmp_yx,1)];
    end
    outFilter.modelMoments = outFilter.resMoments.modelMoments;
    
    if addMerror == 1
        % Adding measurement errors to second moments
        outFilter.resMoments.stdYields   = outFilter.resMoments.stdYields + model.params.stdMr20;
        outFilter.resMoments.stdForEstim = outFilter.resMoments.stdForEstim + sqrt(diag(Rw(1:length(outFilter.resMoments.stdForEstim),1:length(outFilter.resMoments.stdForEstim))));
        outFilter.modelMoments(setupFilter.nyMom-2+1:2*setupFilter.nyMom-2) = ...
            outFilter.modelMoments(setupFilter.nyMom-2+1:2*setupFilter.nyMom-2) ...
            + diag(Rw(1:setupFilter.nyMom,1:setupFilter.nyMom))';
        if setupFilter.SMMwithCSreg == 1 || setupFilter.SMMwithCSreg == 2
            % For CS-moments: Y = y^(k-1)_t+1 + eps^(k-1)_t+1 - y^(k)_t - eps^(k)_t
            %                   = (y^(k-1)_t+1 - y^(k)_t) + eps^(k-1)_t+1 - eps^(k)_t
            %                 X = (y^(k)_t + eps^(k)_t - y^(m)_t - eps^(m)_t)
            %                   = ((y^(k)_t - y^(m)_t) + eps^(k)_t  - eps^(m)_t)
            % Given NID measurement errors we have
            % E(Y*X] = true term - E(eps^(k)_t^2)
            % Var[X] = true term + E(eps^(k)_t^2) + E(eps^(m)_t^2)
            % We may not have individual measurement for all maturities in
            % the CS-regression, so we just use the M-error for the 5-year
            % bond
            numCSmom = length(setupFilter.maturites_CS);
            % For Cov(Y,X)
            outFilter.modelMoments(2*setupFilter.nyMom-2+1:2*setupFilter.nyMom-2+numCSmom) = ...
                outFilter.modelMoments(2*setupFilter.nyMom-2+1:2*setupFilter.nyMom-2+numCSmom) - scaling*model.params.stdMr20^2;
            % For Var(X)
            outFilter.modelMoments(2*setupFilter.nyMom-2+numCSmom+1:2*setupFilter.nyMom-2+2*numCSmom) = ...
                outFilter.modelMoments(2*setupFilter.nyMom-2+numCSmom+1:2*setupFilter.nyMom-2+2*numCSmom) + 2*scaling.^2*params.stdMr20^2;
            % Updating cs loadings to account for measurement errors
            cs_loadings = outFilter.modelMoments(2*setupFilter.nyMom-2+1:2*setupFilter.nyMom-2+numCSmom)...
                ./outFilter.modelMoments(2*setupFilter.nyMom-2+numCSmom+1:2*setupFilter.nyMom-2+2*numCSmom);
            outFilter.resMoments.CS_loadings = [nan cs_loadings];
            
            % Note that if setupFilter.SMMwithCSreg == 2, then no correction
            % for measurement error on risk-adjusted yields can be done
        end
    end

end
   
%% Sign restriction for the real term premium
% setupFilter.setupProj.objectFuncNum = 3; %for REAL bond prices
% modelRealYields      = model;
% modelRealYields      = rmfield(modelRealYields,'thetaBondPrice');
% setupFilterTmp       = setupFilter;
% setupFilterTmp.setupProj.maxMat = 8;   %bond yields only out to 2 years
% modelRealYields      = projectionStepBondPrices(modelRealYields,setupFilterTmp.setupProj);
% setupFilterTmp.maturities = setupFilterTmp.setupProj.maxMat;
% modelRealYields = getTermPremia_loadings(modelRealYields,setupFilterTmp);
% termPremReal = getTermPremia(modelRealYields,outFilter.xhat);
% if mean(termPremReal.TP(end,:),2)<0
%     disp('Negative real term premium')
%     errorMes = 1;
%     return
% end


%% The objective value
T                      = setupFilter.T;
outFilter.momCondition = outFilter.modelMoments-setupFilter.dataMoments;
outFilter.SMMpart      = setupFilter.lambdaSMM*outFilter.momCondition*setupFilter.W*outFilter.momCondition';
outFilter.QMLpart      = -sum(outFilter.LogL(setupFilter.numInitial+1:end))/(T-setupFilter.numInitial); 
outFilter.LogL         = -outFilter.LogL; 
Q                      = outFilter.QMLpart + outFilter.SMMpart;
outFilter.Q            = Q;

% disp('Extra penalty')
% tmp.setupFilter       = setupFilter;
% tmp.model             = model;
% newsDecom             = inflVarianceRatio(tmp,outFilter.xhat);
% Q                     = Q+ nanmean(abs(newsDecom.summaryTable(:,1)-newsDecom.summaryTable(:,2)))*100;
% outFilter.Q           = Q;
% if setupFilter.numSim > setupFilter.burnIn
%     outFilter.outSim              = outSim;
%     outFilter.ySimAsInData        = ySimAsInData;
%     outFilter.ySimAsInDataNoNoise = ySimAsInDataNoNoise;
% end
  
    
end

